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Abstract 

We introduce a Self-affine Asperity Model (SAM) for the seismicity 
that mimics the fault friction by means of two fractional Brownian pro- 
files (fBm) that slide one over the other. An earthquake occurs when 
there is an overlap of the two profiles representing the two fault faces 
and its energy is assumed proportional to the overlap surface. 

The SAM exhibits the Gutenberg-Richter law with an exponent [3 
related to the roughness index of the profiles. Apart from being an- 
alytically treatable, the model exhibits a non-trivial clustering in the 
spatio-temporal distribution of epicenters that strongly resembles the 
experimentally observed one. A generalized and more realistic version 



of the model exhibits the Omori scahng for the distribution of the after- 
shocks. 

The SAM lies in a different perspective with respect to usual models 
for seismicity. In this case, in fact, the critical behaviour is not Self- 
Organized but stems from the fractal geometry of the faults, which, on 
its turn, is supposed to arise as a consequence of geological processes on 
very long time scales with respect to the seismic dynamics. The explicit 
introduction of the fault geometry, as an active element of this complex 
phenomenology, represents the real novelty of our approach. 

PACS NUMBERS: 91.39.Px, 05.45.+j 
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1 Introduction 



Recently, theoretical models have acquired an increasing relevance in the study 
of the seismicity. Their aim is to fill up the gap between the experimental 
knowledge and the theoretical comprehension of the phenomenon. 

One of the most serious problem which geologists have to face is the lack of 
complete catalogues extended over long time periods. This makes difficult to 
improve the general comprehension about earthquakes. By studying theoret- 
ical models, one then tries to focus on some particular ingredients, which are 
supposed to be essentials, and then tries to understand as much as possible of 
the seismic behaviour. In this way one can compare the specific predictions of 
the models with those obtained from the real catalogues. 

Though the dynamics of earthquakes is very complex there are some simple 
basic components which have to be taken into account in a model: 

(a) earthquakes are generated by a very slow discontinuous driving of a fault; 

(b) the occurrence of earthquakes is intermittent, i.e. they occur as abrupt 
rupture events when the fault can no longer sustain the stress; 

(c) there are two separate time scales involved in the process; one is related 
to the stress accumulation while the other, which is orders of magnitude 
smaller, is associated to the duration of the abrupt releases of stress. 

Many forms of scaling invariance appear in seismic phenomena. The most 
impressive feature is the celebrated Gutenberg-Richter law for the magni- 
tude distribution of earthquakes. It states that the probability P{E)dE that 
an earthquake releases an energy in the interval [E,E + dE] scales according 
to a power-law 
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with an exponent j3 of order of unity whose eventual universahty represents 
matter of debate. 

The Omori law for the time correlations of aftershocks (i.e. seismic 
events which happen as a consequence of a main earthquake) is another ex- 
ample of scaling behaviour in the seismic phenomenology and one of the most 
difficult to reproduce in simplified models. 

In the last decades there has been an increasing evidence for the space-time 
clustering |^ of the earthquake epicenters. In particular there are experimen- 
tal evidences suggesting that the epicenter distributions is self-similar both in 
space and in time. 

Unfortunately, the complexity of modelling the motion of a fault system, 
even in rather well controlled situations such as the San Andreas fault in Cal- 
ifornia, is a highly difficult task and it is still controversial what is the correct 
theoretical framework at the very origin of scaling laws. It is thus important 
to individuate models as simple as possible that are able to exhibit the main 
qualitative features of the fault dynamics. Their physical relevance stems from 
the specific predictions on the real seismic activity which might be verified from 
experimental data. 

One of the first attempt in this direction is due to Burridge and Knopoff 
who introduced a stick-slip model of coupled oscillators to mimic the interaction 
of two fault surfaces. In practice, one considers blocks on a rough support 
connected to one another by springs. They are also connected by other springs 
to a driver which moves at a very low constant speed. The blocks stick until the 
spring forces overwhelms the static friction and then one or more blocks slide, 
releasing an 'earthquake' energy proportional to the sum of the displacements. 
In the frame of the inferior plate, if Xi denotes the position of th z-th block. 
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the equations of motion are 

TTliXi = fcc,j(3^i+l ~ ^i) ~ k(. i_i(^Xi — Xi^i) — kp ii^Xi — vt) + Fii^Xi) (2) 

where Fi{xi) represents the friction force which depends on the block velocity 
Xj. In the original model the friction force, zero for zero velocity, increases 
progressively as the velocity increases up to a certain maximum value. This 
model exhibits the Gutenberg-Richter law for the distribution of the energy 
released during an earthquake and it allows for the presence of aftershocks. Up 
to now the original model of Burridge and Knopoff remains the only one able 
to explain the presence of aftershocks without ad hoc modifications. 

A numerical integration of the Newton equations for a one-dimensional 
chain with a large number of homogeneous blocks has been performed by Carl- 
son and Langer [^. Their model differs from that of Burridge and Knopoff 
in the form of the friction force which is supposed to be identical for all the 
blocks, neglecting the inhomogeneities of the crust. It has been shown that the 
model exhibits the Gutenberg-Richter law (see also [^] for the connection 
with the chaotic behaviour of the system). 

More recently it has been suggested that the qualitative aspects of earth- 
quakes (and of Burridge and Knopoff models) could be captured by the so- 
called Sandpile models, which are the paradigm of a large class of models 
showing Self-Organized Critical (SOC) The concept of Self-Organized Crit- 
icality (SOC) has been invoked by Bak, Tang and Wiesenfeld to describe 
the tendency of dynamically driven systems to evolve spontaneously towards a 
critical stationary state with no characteristic time or length scale. An exam- 
ple of this behaviour is provided by Sandpile models: sand is added grain by 
grain in a pile on a c?- dimensional lattice until unstable sand (too large local 
slope of the pile) slides off. In this way the pile reaches a steady state where 
additional sand grains fall off the pile by avalanche events. This steady state 
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is critical since avalanches of any size are observed. According to this pic- 
ture of Self-Organized Criticality, during its whole evolution the earth would 
have reached a marginally stable state in which any small perturbation could 
give rise to relaxation processes, earthquakes in this case, that can be small or 
cover the entire system. In this way the earthquakes would be the equivalent 
of avalanches for Sandpiles models. The main ingredient in this picture would 
be the interplay between the slow dynamics, represented by the stress accu- 
mulation, and the fast dynamics of earthquakes. The latter would modify the 
earth crust which, on its turn, can give rise to earthquakes and so on, with a 
feedback mechanism that would be at the origin of the self-organization. 

There exists a whole generation of SOC models proposed to explain the 
scale-invariant properties of earthquakes ||, [1^ . These type of models suggest 



however that there is no stress accumulation before a big earthquake and the 
exponent of the Gutenberg- Richter law is expected (with the exception |^ that 
we mention hereafter) to be universal. In addition the space-time distribution 
of the epicenters has no clear relation with the experiments where non-trivial 
clustering are present. 

It is worth to recall in this framework the model proposed by Olami, Feder 
and Christensen ||^. Their model maps the two dimensional version of the 
Burr idge-Knop off spring-block model in a cellular automaton and it gives a 
good prediction of the Gutenberg-Richter law with a non-universal value of 
the /3-exponent, which varies with the level of non-conservation of the model 
and could account for the (3 variances observed in nature. 

In order to go beyond the limitations of these models, we have recently 
proposed an alternative approach |11|] where the critical behaviour is not self- 



organized but stems from the fractal geometry of the faults |jT2|, |T4l . In this 



perspective the faults are supposed to be formed as a consequence of geolog- 
ical processes on very long time scales with respect to the seismic dynamics. 
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Looking at the system on the time scale of human records the fault structure 
can be considered assigned and just slightly modified by earthquakes. 

In particular, we have introduced the so-called Self-affine Asperity Model 
(SAM) |TT[ which mimics the fault dynamics by means of the slipping of two 
rough and rigid brownian profiles one over the other. In this scheme an earth- 
quake occurs when there is an intersection between the two profiles. The 
energy released is proportional to the overlap interval. This model, apart from 
being analytically treatable, exhibits some specific features which follow from 
the fractal geometry of the fault. In particular it reproduces the Gutenberg- 
Richter law with an exponent (3 which is non-universal since it depends on the 
roughness of the fault profiles. It predicts the presence of a local stress accu- 
mulation before a large seismic event. Moreover it allows one to analyze and 
investigate the complex phenomenology of the space-time clustering of epicen- 
ters. The model exhibits, in fact, a long-range correlation of the events which 
corresponds to a self-similar distribution of the spatial and temporal epicenter 
sets. In this scheme it is also possible to include the analysis of the origin of 
aftershocks and show that, in a natural generalization of the model, they follow 
the celebrated Omori's law. 

In this paper we describe in detail the SAM. The analytical results are, 
step by step, tested numerically and, whenever possible via the comparison 
with experimental data. 

The outline of the paper is the following. In sect. II we introduce the 
model and we recall some properties of fractional Brownian profiles. Sect. Ill 
is devoted to the discussion of the Gutenberg-Richter law. We show that the 
SAM follows this scaling with an exponent (3 that we relate analytically to the 
roughness of the brownian profile. This allows us to draw some conclusions on 
the non- universality of the exponent (3. In sect. IV we discuss the problem 
of the distribution of epicenters both from the spatial point of view and the 
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temporal one. The SAM exhibits a non-trivial clustering of epicenters which 
reproduces the experimental results and can be analytically explained by ex- 
ploiting the properties of the fractional Brownian profiles. The problem of 
the power spectrum of the temporal sequence of earthquakes is also discussed. 
Sect. V is dedicated to the introduction of a more realistic version of the SAM 
model. This version, which takes into account the local rearrangement of the 
earth crust as a consequence of the earthquakes, exhibits a non-trivial scaling 
in the distribution of the aftershocks, according the Omori's law. Finally in 
sect. VI we draw the conclusions. The paper is completed by two appendices 
on the statistics of the fractional Brownian profiles. 
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2 The model 



Many authors pointed out that natural rock surfaces can be represented by 



fractional brownian surfaces over a wide scale range [T^, [T^ and that also the 
topographic traces of the fault surfaces exhibit scale invariance [O. A fault 
can thus be regarded as a statistically self-affine profile -Fff(t), whose height 
scales as + r) — Fnit)] ~ . In d = 2, such a profile Ffjit) can be 

generated by fractional Brownian motion (fBni) with exponent H, the Hurst 
exponent, and in = 3 by the standard generalization given by brownian reliefs 



T8| , |T9[. The exponent < H < I controls the roughness of the fault where the 
standard brownian profile corresponds to H = 1/2, and a different iable curve 
corresponds to H = 1. Just to give an example let us recall how it is possible 
to generate a brownian profile. In the one- dimensional case one can generate 
L random variables (r.v.) {Xi, ...,Xl} according to the following algorithm: 

1 with prob. p = 1/3 
Xi = { with prob. p = 1/3 
— 1 with prob. p = 1/3 

On a one- dimensional lattice of L sites one can thus define a stochastic function 

S[n) = So + Y.^^ yn<L (3) 

i=l 

where 5*0 is an arbitrary integer number. The (|^) defines, in the limit n ^ oo 
a self-affine profile of fractal dimension D = 1.5. More in general the fractal 
dimension of the profile is well known to he Dp = d — H. For further details 
we refer to the appendix A. 

The explicit introduction of the fault geometry in a model for seismicity 



was already been supposed by Huang and Turcotte ||T2|. They introduced a 
static model where the average of all the seismic events contributing to the 
Gutenberg-Richter law is taken over many uncorrelated realizations of one 
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single fractal profile. The purpose of this letter is to introduce a dynamical 
model, called Self-affine Asperity Model (SAM), that describes the seismic 
activity considering two profiles sliding one over the other instead of only one as 
in |jl2[ . Such a model has the advantage to exhibit strong spatial and temporal 
correlations also between far away seismic events, and allows us to infer some 
specific and new predictions about the relation between the roughness of the 
fault H and the scaling exponent of the Gutenberg-Richter law as well as on 
the spatio-temporal distribution of epicenters. 

Note how this model represents an alternative approach with respect to 
the SOC models. In this case, in fact, one supposes as lacking the interplay 
between the fault structure and the seismic events. The latter is supposed not 
to modify substantially the fault geometry. In this sense one is in a sort of 
limit of infinite rigidity of the Burridge-Knopoff models. 

Operatively, the SAM is defined by the following dynamical rules: 

(i) We consider two profiles, say S'{n) and S"{n), with n = 1, L, on par- 
allel supports of length L at infinite distance. The initial condition is obtained 
by putting them in contact in the point where the height difference is minimal 
so that (see Fig.l) 

S'in) = S"in)+max,^{,_L}{S'ij) - S"ij)} 
n = 1, L 

(ii) The successive evolution is obtained by drifting a profile in a parallel 
way with respect to the other one, at a constant speed v, so that S'{n;t) = 
S'{n-vty, 

(iii) At each time step t, one controls whether there are new contact points 
between the profiles, i.e. whether S'{n;t) — S"{n) < for some x value. An 
intersection represents a single seismic event and starts with the collision of two 
asperities of the profiles. The energy released is assumed to be proportional to 
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the breaking area of the asperities, i.e. the extension of the hypersurfaces, in 
general of dimension {d — 1), involved in the collision of the asperities during 
an earthquake. In the case d = 2 the energy released is given by the sum of 
the lengths of the two segments indicated with A and B in Fig.l; 

(iv) We do not allow the developing of new earthquakes in a region where 
a seismic event is already taking place, i.e., with reference to the Fig.2, we do 
not take into account the earthquakes which eventually take place in the region 
A and B of the two profiles, until A and B have a non-zero overlap. 

Rule III is a consequence of the proportionality between the energy released 
during an earthquake and its seismic momentum Mq, which, on its turn, is 
proportional to the average displacement of a fault during an earthquake. It is 
obviously possible to consider more sophisticated schemes and the work along 
these lines is already in progress. 

With these rules, the motion of the two profiles simulate the slipping of the 
two walls of a single fault. The points of collision are the points of the fault 
where the morphology prevents the free slip: these are the points where there 
is an accumulation of stress and, consequently, a raise of pressure. When the 
local pressure exceeds a certain threshold, it happens a breaking, an earth- 
quake, which allows to relax the stress and redistribute the energy, previously 
accumulated, all around. We assumed that the region between the two sliding 
profiles of the fault is empty or filled by a granular medium, coherently to the 
observation that the fault gauge is a zone of fractured rocks. According to the 



paper by Herrmann et al. one could think to this granular medium as com- 
posed by roller bearings between the two surfaces. The existence of large region 
between the two rough surfaces could then be related to the so-called seismic 
gap, namely an extended area two tectonic plates can creep on each other with- 
out producing either earthquakes or the amount of heat expected from usual 
friction forces. This zone slides and has no influences on the dynamics due to 
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its relatively lower viscosity. 

For sake of simplicity, in this version of the SAM, there is no real breaking 
of the profiles as a consequence of an earthquake and the profiles maintain 
their structures after a crash. Thus are in the opposite perspective than SOC 
models. Since the earthquake dynamics has no effect on the structure of the 
profile. Realistic situations could well correspond to intermediate cases, of 
course. 

It is possible to introduce a more realistic breaking mechanism where there 
is also a modification of the asperity form after an earthquake. We will discuss 
this possibility in sect.V and we will show there how it is possible, in this 
framework, to reproduce the Omori's law. 

It is worth to stress that the SAM exhibits a strong non-locality since a 
collision in a point x, at the time t can trigger, at later time, a subsequent 
event also very far away. One of the main advantage of the SAM consists 
in the possibility of deriving various analytic results using the properties of 
brownian profiles. 

3 The Gutenberg- Richter law and the non-universahty 
of the /3-exponent 

In 1956, Gutenberg and Richter |^ noticed the dependance of earthquakes 
frequency from their magnitude: the greater the magnitude, the smaller the 
frequency. The relation between the frequency and the magnitude of earth- 
quakes is: 




(4) 



where N{M > m) is the number of earthquake with a magnitude greater than 
m while a and h are two empirical parameters. The h value is generally in the 
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range 0.8 < b < 1.4 depending on the earth region considered and the stress 
level of the region itself. 

Relation (^) is the most important statistical representation of seismicity 
and the understanding of the underlying mechanisms is of fundamental im- 
portance for the comprehension of earthquakes and their prevision. Several 
studies have been achieved to understand the origin of the universality of the 
Gutenberg-Richter relation but, despite the simplicity of this relation, there is 
no understanding of the underlying mechanisms. The b value might depend 
on three factor: (1) the geometrical properties of the fault, (2) the physical 
properties of the medium and (3) the stress level of the seismic region. In this 
section we show that, in the framework of the SAM, the b value is essentially 
determined by the fault geometry and in particular by its fractal dimension. 
The magnitude M is not the only indicator of the earthquakes strength; an- 
other quantity used to describe the earthquake intensity is the seismic moment 
defined by the relation: 

Mo= / fiSdA, (5) 

J A 

where /i is the rigidity modulus of the medium under consideration, 5* is the 
slip and A is rupture area. From dimensional analysis it is obvious that the 
energy E released by an earthquake is proportional to its moment. There is an 
empirical relation between the seismic moment (or energy) and the magnitude: 

logio E = cM + d, (6) 

where E is the released energy. From eq. and we easily obtain the 
energy distribution for earthquakes: 

P{E) ~ E~^'\ (7) 

where P{E) is the probability of an earthquake releasing an energy E and 
(3 = b/c. 
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In order to describe the seismic phenomenology a model for the fault slip 
has to verify equation (^: we thus will study the energy distribution for the 
model defined in the previous section (SAM). 

The numerical simulations provide a clear evidence that our model exhibits 
the Gutenberg- Richter law (0), see Fig.3. As we have defined in the previ- 
ous section, the energy released during an earthquake is essentially given by 
the length of the superposition between the fluctuations of the two self-affine 
profiles. Remembering that the difference between two self-affine profiles is 
a self-affine profile itself, we can consider only the profile given by the differ- 
ence between the upper profile and the lower profile: the energy distribution 
will be simply the length distribution of the segment obtained intersecting the 
difference profile with a straight line. 

If we consider a fractal ensemble having a dimension D = d ~ H embedded 
in a (i-dimensional Euclidean space, the intersection between the ensemble and 



an hyperplane of dimension d — 1 will be an ensemble of dimension fl? 

D = {d- H) + (d-l) - d = d- H -1. 

Therefore, the average extension of the hyper-areas given by the intersection 
between a self affine hyper-surface and an hyper-plane will be: 

<a>L~Ad^, (8) 

where the subscript L indicates that we are considering a portion af hyperplane 
of extension A ~ L'^^^. In virtue of the self-affine nature of the considered 
ensemble, the hyper-areas distribution will be: 

d{a) ~ a-'^-\ 

From the last equation we can compute the mean extension of the hyperareas: 

< a a-^da ~ l'-"-')'-^-''^ (9) 
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By comparing and one gets the relation between the exponent (3 of the 
Gutenberg-Richter law in d dimensions and the Hurst exponent which accounts 
for the fractal properties of the faults: 

^^l-j^. (10) 

In the three-dimensional case one has: 

/? = 2-f, (11) 

with/3G [|,1]. 

In order to check the (0) we have performed a numerical experiment in 
d = 2. Fig. 4 reports the results of the /5-value, as a function of the Hurst 
exponent, which are in good agreement with the expected relation (3 = 2 — H . 
The dependence of /3-value on the roughness of the faults could then account 
for the non-universality of the /5-value which would reflect the variability of 
the fractal dimension of the fault profiles around the world. In this perspective 
one could also try to relate the /3-value to the age of a given fault profile. 
By supposing that the effect of the fault slipping and of the earthquakes is a 
smoothing of the profiles, i.e. an increase of iJ, one could guess that the older 
the fault profile, the smaller the /3- value. 

4 Space-time distributions of epicenters 

Let us now try to analyze the problem of the space-temporal clustering of the 
earthquakes epicenters. Many authors |1^, [16|, pointed out that the epicenter 



tend to cover a fractal set with a fractal dimension which is a highly irregular 
function of space and time. One of the most interesting feature is represented 
by the evidence that the spatial distribution of the epicenters along a linear 
seismogenetic structure seems to exhibit self-similar properties. Results of the 
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same kind have been reported for single "transform" faults. This could lead 
to the conclusion that the non-homogeneity of the spatial distribution of the 
epicenters is due to some peculiar phenomenon occurring also in a single linear 
fault and just partly to the fractal distribution of the faults. 

Similar properties are exhibited by the temporal distribution of events whit 
a non-homogenous structure, made of periods of quiescence and bursts of activ- 
ity. Along the same region there could be sub-regions with non-homogeneous 
and also very different behaviours. 

Thanks to the simple dynamics of the SAM model it is possible to study, 
whether analytically or numerically the complex space-time distribution of the 
epicenters. Opcratively the space location of an epicenter is defined in corre- 
spondence of the first point of contact of the two colliding asperities belonging 
to the two profiles. 

As far as spatial distribution of earthquakes is concerned, our simulations 
provide a good evidence of a spatial clustering of epicenters on a set with fractal 
dimension smaller than 1. In particular we obtained a value of the fractal 
dimension dep in the range — 0.8 — 0.9 for if varying in the interval [0.3,0.7] 
and for different lengths of the system between L — 1000 and L — 50000. 

By numerical analysis of the model the values of dgp seem to decrease with 
increasing H and seems to remain nearly constant with respect to variations 
of the system dimensions and of H. These results are not of immediate expla- 
nation; if the fault profiles could slip for an infinite time, in fact, each point 
of the inferior profile could be, theoretically, an epicenter because it, sooner or 
later, would be hit by an asperity of the superior fault profile. In this way we 
would have that limj^oo dgp = limL^oo dep = I, and the set of epicenters thus 
becomes a compact set. This intuitive idea turns out to be correct since the 
observed non-integer fractal dimension is a non-trivial finite-size effect. It is 
possible to show analytically, for H — 0.5, that the fractal dimension dep{L) of 
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the epicenter set in a fault of linear size L is 



dep{L) ~ 1 



7 In In L 



for large L. 



(12) 



InL 



We will sketch here the main lines of the proof referring the reader for further 
details to appendix A. 

First of all it is worth to remind that, according to the definition of the 
SAM, in order to obtain an infinite evolution of the system we necessarily need 
two fault profiles with a length L, which tends to infinity. One has first to create 
the two profiles separately and then to put them in contact. This because the 
average distance between the two profiles tends to increase as L ^ oo. 

In full generality, with reference to Fig. 5, called < S'{n) >= S'q, < S"{n) >= 
Sq and ho = maXjg{i^...^L}{S"(j) — S"{j)}, one has, from eq.(g), < S"{n) >= 
Sq + ho and putting, without loss of generality, S'q = Sq = 0, 



The idea we want to use is that the number of epicenters A'^^, for sufficiently 
large systems, will be proportional to the number of points of the inferior 
profile, N{h,L), with an height h between the minimum value, hmin, of the 
upper fault trace, and the maximum value, hmax, of the lower one, as shown 
in Fig.5: 



where 37]/4 is a constant dependent on the variance of the variables {xi} used 
to generate the profile. By inserting the (^) into the (|1^ one has 



< S"{n) - S'{n) >= h 



(13) 




(14) 



where N{h, L), for big values of h, can be written as 



, 3h,~ 

N{h, L) ~ VLexp 



(15) 




(16) 



17 



Let us find an estimate of tlie hmin and hmax values in tlie limit L ^ oo 
and consider the two faults to be brownian profiles {H =1/2) with length L. 

In our case the variables {Xj} {Yi} which compose the profiles are random 
variables with zero mean and variance = 2/3. So the variables {Xi} = 



3/2{Xi} and {Yi} = y3/2{Yi} will be random variables with zero mean 
and unitary variance. To these variables we can apply the so-called Iterated 
Logarithm Theorem (ILT) It states that, for a partial sum Sk = Y,i=i'^i 
of identically distributed random variables {wi} with < Wi >= and variance 
=< w? >= 1, it holds: 



P lim sup , =1 =1 (17) 

where P{A = a) is the probability for the variable A to have the value a. 
For the H = ^ case we can also write Sinf{k) = J2i=iXi and Ssup{k) = ho + 
J2i=i where {Xi} and {1^} are uniformly distributed variables with zero 
mean, standard deviation cr^ = 2/3 and ho = maxvj(X]fc=i -^i ~ ^i)- 

By using the ILT with profiles built with the normalized variables Xi and 
Yi one obtains: 



limL^ooSup„g[i^i]{S"(n)} = ^vZlnlnL 
limL^ooinfne[i,L]{5"'(n)} = -^VL In InL. 



One has also 



ini {S'{n)} = ho+ inf AS" (n)} (19) 

ne[l,L] ne[l,L] 

By defining the stochastic variables {Zi = Xi — Yi} it will be, by definition, 

n 

ho= sup {Y^Zi}. (20) 

ne[l,L] i=l 

The variables {Zi} have zero mean and variance = 4/3. So we can apply 
the ILT to the variables Zi = ^Zi by getting 

ho = —^VlnlnL. (21) 
V3 
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By comparing (0), (|T9|) and (|2T| ) one easily gets the expression for hmin and 

hmax • 

hmin = ^(V2 - l)VLlnlnL 
hmax = ^vZlnlnL 
and, inserting these expressions in (|T6|) and making a change of variables, 

Ars(L) ~ L / exp-^ dt = LI(L), (23) 

J{'/2-l)VLln InL 

where /(I^) is an integral which tends to zero in the limit L ^ oo. We are 
interested to how this integral goes to zero. 

The " average theorem" for continuous function states that it will be possible 
to find a i = 'j{L) ■ v^2 In InL, with 7(-L) g]v^ — 1, 1[, in such a way that 

where At is the integration interval and 7(L) is the limit value of 'j{L) and we 
have neglected all the terms diverging slower than the logarithm. 
Using the mass-length definition of fractal dimension, 

dep = lim \nNE{L)/lnL (25) 

L^oo 

, we obtain the relation 

^ In In L ^ / In In In L \ , , 

where rj* is the mean value of 7^(L) / r]{L). This implies that, according to what 
we had forewarned, lim^.^oo'^ep = 1, and, thus, that the fractal nature of the 
spatial distribution of epicenters is due to the faults finite size. The asymptotic 
value dep = 1 is reached very slowly at increasing L and it cannot be detected 
but by means of huge simulations. We have checked the validity of (|26|) for 
profiles with a linear size L varying in the range 10^ — 10^. Work is in progress 
to extend our results to the case of a generic roughness index H. 
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Let us now discuss the temporal correlations of earthquakes and in partic- 
ular problem of the 1// noise. A system is said to exhibit 1// noise when its 
power spectrum scales as 

S{f) - (27) 

with a smaller than 2. The interest on 1// noise lies in its ubiquity in nature. 
1/ f noise has been detected in systems as divers as resistors, the hourglass, 
the flow of the rivers or of the cars in a traffic system. Despite much work has 
been devoted to this topic it is still lacking a general theory that explains the 
widespread occurrence of 1// noise. 

The fact that the power spectrum is connected to the autocorrelation func- 
tion by the Wiener and Khintchin theorem leads some authors to the idea that 
the presence of the 1// noise indicates the presence of self-similarity in the 
distribution of correlation times. 

The autocorrelation function is usually defined as 

where E{t) in our case represents the energy released by an earthquake occurred 
at the time t and the averages are taken over the distribution of times to- If the 
energy presents a power-law distribution with an exponent greater than —2, 
as in our case, the average < E(t) > will depend on its maximum value and 
then on the system dimension. We would have, in this way, a non consistent 
procedure to calculate the autocorrelation function. In order to overcome his 
difficulty one can use an alternative definition of the autocorrelation function 
which is independent of the scale of the system. If we define it as 

C{t) =< E{t + to)E{to) > (29) 



it is possible to show |j22| that the power spectrum S{f) is linked to the Fourier 
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transformation of C(t), < l-E/P >, by the relation 

KlEfl"^ >= S{f) + 1/N (30) 

where N represents the dimension of the system. We have also 

S{f)=j:<\Ef'\'>\W{f-f)\'. (31) 
/' 

where the function W{f) takes into account the finite dimension of the system 
and tends to the delta function 5(0) for an infinite system. If is big enough 
one has: 

S{f) ^< \Ef\' > (32) 

and one can study the power spectrum simply analyzing the Fourier transform 
of the autocorrelation function (P5|). 

In our numerical simulation we have studied this function and the results, 
shown in Fig. 6, gave 

S{f) ~ (33) 

with a ~ 1.2. This means that our model exhibits 1/f noise, i.e. it does not 
exist a maximum autocorrelation time and a seismic event may be influenced 
by another one very distant in time. 

5 The Generalized SAM model 

Up to now we have studied a version of the SAM model corresponding to the 
limit of infinite rigidity of the faults. The fault profiles are not modified by the 
seismic activity and one studies the statistics of earthquakes in the hypothesis 
that there is a complete time-scale separation between the seismic activity and 
the rearrangement of the earth crust. The latter would develop in very long 
times with respect to the scale of human records and this would justify the 
assumption. 
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We have shown how this model exhibits a good interpretation of the seismic 
phenomenology in a global sense: Gutenberg- Richter law, epicenter clustering. 
What is lacking is the description of what happens locally, i.e. as a consequence 
of a single event, both from the temporal point of view or from the spatial one. 
In particular it is not possible to obtain in such a scheme the Omori's law for 
the distribution of aftershocks. These events are related to the situation in 
the neighborhood of the main shock epicenter after the occurrence of the main 
shock. They are ruled by the following empiric relation 



where N{t) indicates the number of earthquakes occurred at the time t after 
the main shock, c is a constant and a is an exponent whose value ranges in the 
interval [1.0 -r- 1.4]. For enough long times t one usually supposes t » c and 
the functional form of N{t) is given by a pure power-law N{t) ~ 

In this section we improve the model in order to include the rearrangement 
of the earth crust as a consequence of the occurrence of an earthquake. With 
this modification it is possible to describe the local phenomenology of seismicity 
and in particular to reproduce the Omori's law. 

The model is modified by considering the asperity breaking in the collisions. 
When two asperities collide a fracturing process starts in the smallest asperity 
(that one with the smallest section at the level of the epicenter). The fracture 
propagates inside the fault until it crosses again the fault profile. At this 
point the fracture stops and the resulting configuration represents the new 
fault profile in the region interested by the earthquake. The magnitude of 
the earthquake is assumed to be proportional to the linear extension of the 
fracture. 

In Fig. 7 it is shown an example of fracturing process during an earthquake. 




1 



(34) 




22 



The shadowed region is removed from the fault profile. The statistical proper- 
ties of the fracture are supposed to be identical to the ones of the entire fault 
profile. This means that one has to consider a self-affine profile with the same 
Hurst exponent of the original fault. 

In our simulations, we considered, for sake of simplicity, the case of a Brow- 
nian profile with H = 0.5. Let us note that an earthquake in a certain point 
can trigger several other earthquakes, with smaller magnitude, which occur in 
the same region or in a very close region. In order to investigate the statis- 
tics of the aftershocks we identified all the aftershock occurring after a certain 
mainshock in the rupture region. A mainshock is defined as an earthquake 
above of a certain magnitude (in our simulations an earthquakes involving at 
least 100 sites). Starting from this event one counts, as a function of the time 
elapsed from the mainshock, the number of earthquakes, with a magnitude 
smaller than that of the mainshock, occurring in the same region. One stops 
the counting when an earthquake with a magnitude greater or equal to the 
mainshock occurs. 

We have studied the behaviour of the cumulative distribution of aftershocks, 
i.e. the number Ncum.it < T) of earthquakes occurring before T time steps after 
the mainshock. By averaging over many realizations (of the order of 10^) we 
have obtained the curve reported in Fig. 8 that exhibits the Omori scaling law 



(|3^ . For values of t large enough one has the power law: 

Ncum{t < T) ~ T^-". (35) 

with the exponent a ^ 0.37. The numerical value of the exponent a is not in 
good agreement with real values. However, we have just considered the case 
of a one- dimensional profile embedded in a two-dimensional space and the 
model considers only one isolated faults, thus neglecting the effects of inter- 
action among different faults. It would be interesting to study what happens 
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considering the case of a two-dimensional surface too. 

This generahzed SAM recalls the work of Herrmann et al. about the 
space-filling bearing. The analogy lies in the fact that one could think the 
interspace between the two fault planes as filled by a granular medium which 
is also composed by the broken asperities of the fault. The link is made closer 
by the fact that in our case the distribution of areas of the asperities broken 
follows a power-law 

P{Aasp) ~ A-^^ (36) 

with an exponent a which could be related analytically to the roughness ex- 
ponent by the relation 

Relation ( p6D is obtained by supposing that the area of the broken asperi- 
ties scales with its linear extension / as Aasp ~ l^^^ by a standard variable 
change. Work in this direction is still in progress and we refer the reader to a 
forthcoming paper p3| . 



It is obviously possible to consider more realistic generalizations of the 
breaking mechanism, in which the application of the pression in a certain point 
causes the breaking in a different point, mimicking, in this way, the effect 
of the stress redistribution in the medium. This situation is, on its turn, a 
simplification with respect to the ideal case in which one has to calculate, at 
each time step, the new stress field in the whole medium as a consequence of 
the changed pression conditions. 



6 Conclusions and perspectives 

In summary, we have proposed a model of earthquakes where the critical behav- 
ior is generated by a pre-existent fractal geometry of the fault. The statistics 
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of earthquakes is thus related to the roughness of the fault via the scaling re- 
lation (2) between critical indices. This result suggests that the younger the 
fault system, the larger the b exponent is, since one expects that the roughness 
of a fault decreases in geological times. Note that in this case, the exponent 
b is non-universal. Another major result is that the fractal distribution of the 
epicenters could be a finite size effect very difficult to be detected from data 
analysis. In our case our results provides a possible explanation for the highly 
irregular and non random distribution of epicenters that is observed experi- 
mentally. Least but not last, the accumulation of pressure is at the very origin 
of large seismic events in the SAM. The presence of such an effect could be 
tested also in real situations e.g. by piezo-electric measurements. 

Moreover we introduced a generalization of the SAM which includes the 
effect of the breaking of the asperities in contact during an earthquake. This 
makes the model much more realistic and allows for the interplay between 
earthquakes and structural properties of the faults. This version of the model 
exhibits a non-trivial distribution of aftershocks which follows the Omori's law. 
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Appendix A: Statistics of fractional Brownian motions 

In this appendix we review the main properties of the so-called fractional 
Brownian motions (fBm for short) which represent a generalization of the Brow- 
nian motion 0, |2|]. 

A fBm Fnit) is defined as a monodrome function of one variable t, such 
that its increment AFniAt) = Fnit + At) — Fnit) has a gaussian distribution 
with variance 

=< AFff^iAt) >~ At^^, (A. 1) 

where the brackets indicate the average over many realizations of Fnit). The 
parameter H is the so-called Hurst exponent and takes values between and 
1. The main properties of those functions can be summarized as follows: 

1) they are stationary, i.e. the average square increment depends only on 
the increment of the argument t and all the values of this argument are 
statistically equivalent; 

2) they are continuous functions but nowhere differentiable; 

3) they are self-affine curves, i.e. if the time scale is rescaled by a factor r, 
the corresponding increment AF^it) is rescaled by a factor r^: 

< AFu^rAt) >~ r^" < AFu^At) > . (A. 2) 

THe fBm are self-affine curves which present a box-covering dimension equal 
to Dp = d — H. Let us consider, for sake of simplicity, the case d = 2 and 
suppose that Fnit) is defined in a time interval At = 1 with a vertical extension 
AFfjlt) = 1. If one rescales the time by a factor r < 1, then, by virtue of the 
self-affinity, F^it) will be rescaled by a factor . Thus, in order to cover a 
section of curve extending in the interval At = r one needs AFn/At = r^^^ 
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boxes of linear dimensionr and for the entire profile one will need ^ jr boxes. 
So recalling the definition of box-covering dimension 

D, = limi^^ (A. 3) 

r-»o logl/r 

one has 

loe r^~^ 

Dn = lim -2 — — (A. 4) 

r-o logl/r ^ ^ 

In the general d-dimensional case one can define the brownian hyper-surface 
as a function of n = d — 1 variables Xj i = 1, n such that 

< AFH^Ar) >~ Ar^", (A. 5) 

with Ar^ = AXi^ + ... + AX^^. 

The box-covering fractal dimension is then defined as 

Dp^n + l- H ^d- H. (A. 6) 

A last word about the intersection of a fBm with a line parallel to the tem- 
poral axis (fractal dimension Di) and lying in the same plane of the brownian 
profile (fractal dimension D2). In this case, by using the law of additivity of 
the CO dimension, 

Do^Di + D2- d, (A. 7) 

where Dq is the fractal dimension of the intersection set, the zeroset. In our 
case one has Di — d — 1 and D2 — d — H and then 

Do = d-l-H. (A. 8) 

The set of zeroes of a Brownian profile in c? = 2 with a generic vale of the Hurst 
exponent H is then a set of points whose fractal dimension is Do = 1 — -f^ 
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Appendix B 

In this appendix we calculate the number of points that, in a Brownian 
profile, lie at a certain height h. from the general properties of the brownian 
profile one knows that if /i = 0, this number is proportional to y/L where L is 
the length of the profile. Moreover, as a consequence of the spatial homogeneity 
of the random walk one has 

P{Sn+^ = 0|5„ = 0) = P{Sn+m = h\Sn = h) (B. 1) 

where P{a\h) is the conditional probability that, given a certain event 6, the 
event a occurs. The number of points at the height h will be proportional to 
\/ L — t where t is the first passage time at the height h. The first passage time 
distribution for an height h is known to be 

One then has that the number of points at the height h is 

N{h) ~ ["^ VL^fhit)dt. (B. 3) 

Jo 

Eq.( p. 3]) is a very complicate expression and we limit ourselves to consider 
what happens just in the range h > a/I. For the average theorem it will exist 
a value t* such that 

N{h) ~ VL^h{e)L. (B. 4) 

We are interested in the case of /i ~ \fL and we can then suppose t* ^ riL. 
One obtains 

N{h) ~ Vlexp-^ (B. 5) 

that we used in sect. 4. 
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FIGURES 




Figure 1: Fault planes realized by two Brownian profiles put in contact in one 
point. 
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Figure 2: Sketch for the definition of the energy released during an earthquake. 
It is assumed proportional to the breaking area (the {d — l)-dimensional sets 
A and B) between the two asperities: E oc A + B. 
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Figure 3: Probability density of the earthquakes releasing an energy E vs. E 
for roughness index H = 1/2. 
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Figure 4: Exponent /? + 1 vs. the Hurst exponent H for the SAM m d = 2. 
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Figure 5: Scheme illustrating the region of heights [hmin, hmax] in which it is 
possible the occurrence of collisions between asperities, ho is a function of the 
length L of the profiles and indicates their average distance. 
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Figure 6: Power spectrum (solid line) for the temporal sequence of earthquakes 
in the SAM model. It shows a 1// behaviour with an exponent a ~ 1.2 
corresponding to the slope of the dashed line. 
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Figure 7: Scheme illustrating the mechanism for the breaking of the asperities 
in the generalized SAM model. When there is a coUision between two asperities 
the weaker is broken. The shadowed region defines the broken area and the 
new profile after the collision. 
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Figure 8: Cumulative distribution for the aftershocks: N{t < T) is the number 
of aftershocks, events causally connected to the mainshock, occurred up to the 
time T, elapsed from the mainshock time. 
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